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We investigate a general scheme for generating, either dynamically or in the steady state, continuous variable 
entanglement between two mechanical resonators with different frequencies. We employ an optomechanical 
system in which a single optical cavity mode driven by a suitably chosen two-tone field is coupled to the two 
resonators. Significantly large mechanical entanglement can be achieved, which is extremely robust with respect 
to temperature. 
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I. INTRODUCTION 

Entanglement is the distinguishing feature of quantum me¬ 
chanics and is the physical phenomenon according to which 
only the properties of the entire system have precise values, 
while the physical properties of a subsystem can be assigned 
only in reference to those of the other ones. It is now inten¬ 
sively studied because it corresponds to peculiar nonlocal cor¬ 
relations which allows performing communication and com¬ 
putation tasks with an efficiency which is not achievable clas¬ 
sically m. 

Furthermore, for a deeper understanding of the boundary 
between the classical and quantum world, it is important to 
investigate up to which macroscopic scale one can observe 
quantum behavior, and in particular under which conditions 
entanglement between macroscopic objects, each containing 
a large number of the constituents, can arise. Entangle¬ 
ment between two atomic ensembles has been successfully 
demonstrated in Ref. EL while entanglement between two 
Josephson-junction qubits has been detected in Refs. 012. 
More recently, macroscopic entanglement has been demon¬ 
strated in electro-mechanical systems 0: continuous vari¬ 
able (CV) entanglement, similar to that considered by Ein¬ 
stein Podolski and Rosen (EPR) 0, has been generated and 
detected between the position and momentum of a vibrational 
mode of a 15 /rm-diameter A1 membrane, and the quadratures 
of a microwave cavity field, following the theory proposal of 
Ref. Q. 

Entanglement between two mechanical resonators (MRs) 
has been instead demonstrated only at the microscopic level, 
in the case of two trapped ions 0, and between two single¬ 
phonon excitations in nano-diamonds 0. The realization of 
this kind of entanglement at the more macroscopic level of mi¬ 
cromechanical resonators would be extremely important both 
for practical and fundamental reasons. In fact, on the one 
hand, entangled MRs at distant sites could represent an im¬ 
portant building block for the implementation of quantum net¬ 
works for long-distance routing of quantum information US); 
on the other hand, these nonclassical states represent an ideal 
playground for investigating and comparing decoherence the¬ 
ories and modifications of quantum mechanics at the macro¬ 
scopic level | [TTUr3l . 


Many different schemes have been proposed in the litera¬ 
ture for entangling two MRs, especially exploiting optome¬ 
chanical and electromechanical devices M EG), in which the 
two MRs simultaneously interact with one or more electro¬ 
magnetic cavity fields. Refs. EMU) considered the steady 
state of different systems of driven cavities: Ref. fThfl focused 
on two mirrors of a ring cavity, while Ref. 03 assumed to 
drive two independent linear cavities with two-mode squeezed 
light transferring its entanglement to the cavity end-mirrors. 
Ref. ESI instead considered a double-cavity scheme in which 
one cavity couples to the relative motion of two MRs, and the 
second cavity to their center-of-mass; when the system is ap¬ 
propriately driven by squeezed light, such squeezing is trans¬ 
ferred to the two MRs which are then prepared in a station¬ 
ary EPR-like state. Actually, steady-state entanglement can 
be achieved, even if at a smaller value, also without squeezed 
driving, either between two movable mirrors in a Fabry-Perot 
cavity ED, between two mechanical modes of a single mov¬ 
able mirror 1201 , or in the case of two semi-transparent mem¬ 
branes interacting with two driven cavity modes H21 1 . 

A different approach for generating entangled MRs exploits 
conditional measurements on light modes entangled or corre¬ 
lated with mechanical degrees of freedom 12214271 . In this 
case, entanglement is generated at the measurement and it 
has a finite lifetime which may be severely limited by the in¬ 
teraction of the MRs with their reservoirs. A similar strat¬ 
egy has been provided to enhance the entanglement of two 
MRs ll28l . More recent proposal applied reservoir engineer¬ 
ing ideas l29ll33l to optomechanical scenarios, by exploiting 
suitable multi-frequency drivings and optical architectures in 
order to achieve more robust generation of steady state entan¬ 
glement between two MRs 044440) . eventually profiting from 
mechanical nonlinearities and/or parametric driving 14T1142 I. 

In the present paper we propose a novel optomechani¬ 
cal/electromechanical scheme for the generation of remark¬ 
ably large CV entanglement between two MRs with differ¬ 
ent frequencies, which is also extremely robust with respect 
to thermal noise. The scheme is particularly simple, involv¬ 
ing only a single, bichromatically-driven, optical cavity mode, 
and optimally works in a rotating wave approximation (RWA) 
regime where counter-rotating, non-resonant, terms associ¬ 
ated with the bichromatic driving are negligible. The scheme 
shares some analogies with the reservoir-engineering schemes 
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of Refs. Il34l[36ll38ll40l . but it may be used to generate robust 
entanglement also in a pulsed regime, in the special case of 
equal effective couplings at the two sidebands, where the sys¬ 
tem becomes analogous to the Sprensen-Mplmer scheme for 
entangling trapped ions in a thermal environment li43l . This 
latter scheme has been already considered in an optomechan¬ 
ical scenario by Kuzyk et al. 1441 for entangling dynamically 
two optical modes via their common interaction with a single 
MR. 

The paper is organized as follows. In Section II we de¬ 
rive the effective quantum Langevin equations (QLE) describ¬ 
ing the dynamics of the system in the RWA. In Section III 
we solve the dynamics in terms of the mechanical Bogoli- 
ubov modes of the system 134] [36i l45l . derive the steady 
state of the system in the stable case, and provide simple an¬ 
alytical expressions for the achievable mechanical entangle¬ 
ment, showing its remarkable robustness with respect to tem¬ 
perature. In Section IV we instead consider the special case 
of equal couplings, when the system can be mapped to the 
Sprensen-Mplmer scheme |43l , in which mechanical entan¬ 
glement is generated only dynamically and slowly decays to 
zero at long times. In Section V we solve and discuss the ex¬ 
act dynamics of the system in order to establish the conditions 
under which the RWA does not seriously affect the robust gen¬ 
eration of large mechanical entanglement. In Section VI we 
discuss the experimental detection of such entanglement and 
present some concluding remarks. In the Appendices we pro¬ 
vide some detail on the dynamical evolution of the system, 
and present a careful derivation of the linearized QLE in the 
RWA regime. 


II. SYSTEM HAMILTONIAN AND DERIVATION OF THE 
EFFECTIVE LANGEVIN EQUATIONS 

As shown in Fig. |T] we consider an optical cavity mode 
with resonance frequency at c and annihilation operator a in¬ 
teracting via the usual optomechanical interaction with two 
different MRs, with frequencies a>] and wi and annihilation 
operators b\ and b 2 respectively. The cavity mode is bichro- 
matically driven at the two frequencies cjq + o>\ and ojq - a> 2 , 
with the reference frequency cjq detuned from the cavity res¬ 
onance by a quantity A 0 = oj, - a> 0 . If we describe the cavity 
field in a reference frame rotating at the frequency ojo , then 
the system Hamiltonian is given by 

H = hco\b\b\ + ha> 2 b\b 2 + hAoa^a 

+h [g! (bi + £{) + g 2 (b 2 + £j)] 

+h [(Lie-'" 1 ' + E 2 e iW2 ') + H.C.] . (1) 

This means that the cavity mode is simultaneously driven on 
the blue sideband associated with the MR with annihilation 
operator b\, and on the red sideband associated with the MR 
with b 2 - The nonzero detuning A (i makes the present scheme 
different from the one studied in the supplementary material 
of Ref. ||34| which restricts to the resonant case A 0 = 0. Our 
model is instead related to the scheme proposed by Kuzyk 


et al. Ii44l for entangling dynamically two optical modes via 
their common interaction with a single MR: here we will dy¬ 
namically entangle two MRs via their common interaction 
with an optical mode. 

The system dynamics can be efficiently studied by lineariz¬ 
ing the optomechanical interaction in the limit of large driv¬ 
ing field. In this case the average fields for both cavity, a(t), 
and mechanical degrees of freedom, flj(t), are large, and one 
can simplify the interaction Hamiltonian at lowest order in the 
field fluctuations 


5a(t) — a(t ) - a(t), 

Sbj(t) = bj(t)-/3j(t). (2) 


Differently from the typical optomechanical settings in which 
the steady state average fields are time-independent, here the 
bichromatic driving induces a time-dependent, periodic steady 
state average field which, in turn, implies time-dependent 
effective coupling strengths for the linearized dynamics of 
the fluctuations. As originally discussed in l46l . and de¬ 
tailed in Appendix |B| approximated dynamical equations for 
the fluctuation operators bait) and 6bj(t) can be derived, in 
the interaction picture with respect to the Hamiltonian //(> = 
h[a>\b\b\ + oj 2 b\b 2 ), by neglecting the non-resonant/time- 
dependent components of the effective linearized interactions. 
It is possible to prove that this approach is justified when (see 
Eq. pTT)) 


Sj 


, K U)j,\(jJ\ - U>2\ ■ 


(3) 


The corresponding QLE including thermal noise and dissipa¬ 
tion at rates k and y ; for the cavity and the mechanical mode 
j e 1,2 respectively, are 


da - — (k + /A) 5a - iG\5b' x - iG 2 5b 2 + V2 ~Ka m (4) 
5b\ = --y 6b\ - iG\6a' i + VyTfej n > (5) 

6b 2 = -jSb 2 - iG\5a + yjyffi, (6) 

where 


G i 
G 2 


giE\ 


u)\ - A + Ik 

g2 (-2 

a> 2 + A — iic' 


(7) 


are the (generally complex) linear optomechanical cou¬ 
plings, and a"' and b'J are standard input noise operators 
with zero mean, whose only nonzero correlation functions 
are (a in (t)a in (t , ) t ) = 6(t - f), (b'J(t) b'"(t'y) = (Fij + 
1 )6(t - t’) and (b'"(t) f b'”(t r )j = hj5(t - t'), where Fij = 

[exp (hu)j/k B T^ - lj is the mean thermal phonon number 
of the y-th MR, which we assume to stay at the same en¬ 
vironmental temperature T. Moreover, we note that here 
the new cavity detuning A includes the time-independent fre¬ 
quency shift induced by the optomechanical interaction, pro¬ 
portional to the DC component of the average mechanical os¬ 
cillation amplitude /3j(t), that we here denote with (see 
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FIG. 1: Sketch of the proposed entanglement generation and detection scheme (a), and of the various pump and probe laser frequencies (b). 
The cavity mode is bichromatically driven at the two frequencies a> 0 + at, and u> 0 - u> 2 - Large and robust entanglement of the two mechanical 
resonators can be generated either dynamically or in the steady state. Two weak probe fields with detuning A'' = ut C j - co 1 ’ = coj, j = 1, 2, are 
then sent into the cavity. By homodyning the probe mode outputs, the mechanical quadratures (xj, Pj ) are therefore measured, which allows 
one to construct the correlation matrix of the quadratures from which entanglement can be derived in a straightforward way. 


Appendix [B). Specifically 

A = A 0 +2]T ^Ref/lf] • (8) 

7 = 1,2 


We will see that the dynamics described by these equations 
allows to generate large and robust entanglement between the 
two MRs, either in the steady state or, in a particular parameter 
regime, during the time evolution with a flat-top pulse driving. 
We first notice that the system is stable when all the eigenval¬ 
ues associated with the linearized dynamics of Eqs. 0-0 
have negative real parts. The stability condition is quite in¬ 
volved in the general case, but it assumes a particularly simple 
form in the case of equal mechanical dampings, y, = yi = 7- 
In such a case, the system is stable if and only if 


|G 2 | 2 >|Gil 2 - 


kj 

2~ 


4A 2 

(r + 2 k) 2 


(9) 


This stability condition reduces to the one derived in the sup¬ 
plementary material of Ref. Il34l in the case A = 0. We see 
that a nonzero detuning generally helps in keeping the system 
stable. 


III. DARK AND BRIGHT BOGOLIUBOV MODES 

The coherent dynamics corresponding to the Eqs. 0-0’ 
is described by the effective linearized Hamiltonian 

H e fi = HAda^da + h(G\6b\ + G^Sb^da^ 

+h(G\6b x +G* 2 6bl)6a. (10) 


We can always adjust the phase reference of each MR (which 
will be determined by a local oscillator which must be used to 
measure the mechanical quadratures for verifying entangle¬ 
ment) so that we can take both G\ and Go real. 

Eq. {lO} naturally suggests to introduce two effective me¬ 
chanical modes allowing to simplify the system dynamics. We 
assume for the moment Gi > G\, which is a sufficient condi¬ 
tion for stability (see Eq. ([9])), and define 

G 1 5b { +G l 6bX „ 

Pi = -—--= 6 b\ cosh r + 6 b ' 2 sinhr, (11) 

„ G2&bn + Gi6b\ „ A + 

fio = - : -= 6 t >2 cosh r + db. sinhr, (12) 

Q 

where 

&= ^ tanhr=pi. (13) 


Eqs. (11 1 -( 12 1 define a Bogoliubov unitary transformation of 
the mechanical mode operators, which can also be written as 


0 l2 - e - r ( S *'l s *i- s & lS h)db l 2 e r (^lsbl-Sbidb 2 ) 

- S(r)6b l2 S(-r), (14) 

with S ( r ) the two-mode squeezing operator. The Bogoliubov 
mode Pi describes the “mechanical dark mode”, which does 
not appear in H e s, i.e., is decoupled from the cavity mode 
and therefore is a constant of motion in the absence of damp¬ 
ing, while P 2 is the “bright” mode interacting with the cavity 
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mode. This is equivalent to say that the dark mode j3\ is the 
normal mode of the Hamiltonian dynamics with eigenvalue 
equal to zero. The other two normal modes of the system will 
be linear combinations of $2 and da. The Bogoliubov mode 
description has been already employed in cavity optomechan- 
ics, associated to two optical modes in Refs. mm, and to 
two mechanical modes in Refs. |36il38fl (see Appendix A for 
a derivation of the normal modes of the system and a study of 
its Hamiltonian dynamics). 


with the number of excitation of the cooled bright mode and 
the correlations between the two Bogoliubov modes respec¬ 
tively given by 


n C 2 °\r) = nf(r) 



(1 ~ e) C 
1 + 6 2 + C 


(23) 


and 


nifj(r) = m(r ) 


2(1 + id) 

2(1 + id) + (1 - e) C ’ 


(24) 


A. Stationary entanglement for different couplings 

For a realistic description of the system dynamics we must 
include cavity decay and mechanical dissipation. It is conve¬ 
nient to rewrite the QLE in terms of the Bogoliubov modes, 
which in the case when 71 - 72 = 7 assume the simple form 


da - ~(k + /A) da - iQfii + V2 Ka m (15) 

h VryS'r, (16) 

k = -\h- i&Sd + (17) 

where jf , j — 1,2, are two correlated thermal noise operators 
whose only nonzero correlation functions are 

= \nf(r)+\\d(t-t"), (18) 


(pyfpy 1 )) - nfir)d(t-t'), 

(p'"(t)'$(/)) = (Jf(t) f = m(r)5(t - /), 

with the effective mean thermal phonon numbers 

hf(r) = n\ cosh 2 r + (IF + 1) sinh 2 r, (19) 

hf(r) = ni cosh 2 r + (hi + 1) sinh 2 r, (20) 

and the inter-mode correlation 

fh(r ) = coshr sinh r(Ji\ + m + 1) . (21) 

If 71 + 72 a dissipative coupling term between the two Bogoli¬ 
ubov modes appears, which however does not have relevant 
effects because it is proportional to |yi - 72 I which is typically 
very small with respect to all other damping rates. 

The dynamics associated with Eqs. (fT5j)-([T7j) is simple: the 
bright mechanical mode jh_ is cooled by the cavity, while the 
correlated reservoir create finite correlations between dark and 
bright modes. In particular, the matrix of correlation for the 
vector of operators fi = {p\,Pi,P\\P 2 }\ whose elements are 
{C/?}^ = (\P)j {/?}<;) is given, at the steady state, by 


C. S = 


' 0 mp nf + 1 O' 
rhp 0 0 «“ o1 + 1 

0 0 trip 

K 0 n“ ,oi m* 0 


( 22 ) 


where 


2 Q 1 

(25) 


jk 


7 

(26) 

7 + 2k’ 

2A 

(27) 

7 + 2k 


and C can be seen as an effective collective optomechani¬ 
cal cooperativity. The steady state correlation matrix can be 
expressed in terms of the original modes b\ and by invert¬ 
ing the Bogoliubov transformation introduced in Eqs. HU> and 
m The result is 


C b = UCpU T 


0 fhb hbi + 1 O' 
iiib 0 0 hb 2 + 1 

hbi 0 0 fhb 

.. 0 h b 2 mb 0 


with 


U = 


cosh r 0 

0 cosh r 
0 - sinh r 

- sinh r 0 


0 - sinh r ' 

- sinh r 0 
cosh r 0 
0 coshr , 


and where now 


hbi = «i ff + sinh 2 r (l + 

-2 cosh r sinh r Re , 
h b 2 - if™ 01 + sinh 2 r (l + nf + « 2 ° o1 ) 

-2 cosh r sinh r Re (hifj , 
ifib = cosh 2 r hip + sinh 2 r hip 

- cosh r sinh r (l + hf + /l™ 01 ) . 


(28) 


(29) 


(30) 


The entanglement between modes b\ and hi, measured by 
means of the logarithmic negativity B71 j48l . can be easily 
expressed in terms of these matrix elements as ll49ll 


En = max [0, - In v], 

v = 1 + h b \ + h b2 - -y4 \m b \ 2 + (h b 1 - h b i ) 2 . (31) 

When the collective cooperativity C is sufficiently large, i.e., 
C » hi(r ), then fhp(r) is negligible (see Eq. (|24|i). This is the 
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working regime in which we are particularly interested, be¬ 
cause in this case, the second Bogoliubov mode can be cooled 
close to its ground state («“ ol (r) <K n| ff (r)), corresponding 
to an entangled state for the original mechanical modes. In 
this case the steady state correlation matrix for the Bogoli¬ 
ubov modes, in Eq. \22\ , reduces to the correlation matrix of 
a state given by the product of two thermal states with occu¬ 
pancies and h^ K>[ (r) respectively. For the two MR of 

interest, associated with the operator b\ and b 2 , such a state is 
just a two-mode squeezed thermal state EDI 


P 1.2 - ^ ( r )P«f (r),th ® Pnf ol (r),th^( -r )> 


where S ( r ) is given in Eq. (141, and 


(32) 


Pn, th 


=2 

n =0 


(1 + ri) r 


^\n)(n\ 


(33) 


is the density matrix of the thermal equilibrium state of a res¬ 
onator with occupancy h. Such a state is entangled for suffi¬ 
ciently large r and not too large mean thermal excitation num¬ 
ber. 

This prediction of large stationary entanglement is con¬ 
firmed in Fig.[2j where we plot the time evolution of the entan¬ 
glement between the two MRs, quantified in terms of the loga¬ 
rithmic negativity En, obtained from the solution of Eqs. (|4j- 
Figure 2 refers to an experimentally achievable set of 
parameters, y = 10 s' 1 , k = 10 5 s' 1 , G 2 = 10 5 s -1 , A = 10 3 
s' 1 , and to different values of mean thermal phonon numbers 
h\, n 2 , and of the ratio G \ / G 2 ■ We see that remarkable values 
of E n are achieved at low temperatures, and that stationary 
mechanical entanglement is quite robust with respect to tem¬ 
perature because one has an appreciable value of E N — 0.32 
even for h\ = 2000, T 12 = 1000. The time to reach the steady 
state is essentially given by the inverse of the cooling rate of 
the bright Bogoliubov mode, which is approximately given by 
^ ^ (/c 2 + A 2 )/(^ 2 /c) Csee Eqs. 

Eq. ( [32] ) suggests that one could achieve large stationary en¬ 
tanglement between the two MRs by taking a large two-mode 
squeezing parameter r, and a large collective cooperativity 
C » 1 in order to significantly cool the bright Bogoliubov 
mode. However the corresponding optimization of the system 
parameters, and especially of the two couplings G 1 and G 2 , 
is far from being trivial. In fact, r increases when G\ —> G 2 , 
which however implies, at a fixed value of G 2 , a decreasing 
value of Q and therefore of C (see Eq. ( fT3) > and Eq. ([25}); 
moreover increasing r has also the unwanted effect of increas¬ 
ing ihp(r) that is the correlations between the two Bogoliubov 
modes (see Eqs. <[21} and <[24])). 

However, a judicious choice of parameters is possible, al¬ 
lowing to get very large stationary mechanical entanglement, 
even in the presence of non-negligible values of the thermal 
occupancies h\ and iij- At a given value of Gi, this is obtained 
by taking a sufficiently large value of the associated single¬ 
mode cooperativity, C\ = 2 G^/kj » 1, and correspondingly 
optimizing the value of G 2 , i.e., of r. In fact, the logarithmic 
negativity associated with the stationary state of Eq. ([32]) can 


be evaluated in terms of the parameter 


Im^O 


[n+(r) + 1] (cosh 2 r + sinh 2 r) (34) 

- y]rt-(r ) 2 + 4 [Fi+(r) + l] 2 sinh 2 r cosh 2 r , 


where n ± (r ) = «j ff (r) ± n“ ol (r), and «“ Kjl (r) can be explicitly 
rewritten in terms of the cooperativity C\ as 


«“ ol (r) = [«2 cosh 2 r + (n\ + 1) sinh 2 rj 
x [ 1 (1 -e)Ci 

sinh 2 r (1 + 6 2 ) + C\ 


(35) 


The dependence of /2 V versus r, for given values of Ci, hi and 
« 2 , shows a maximum and then decays to zero for large r (see 
Fig.[3]which refers to C\ = 2 x 10 4 and rt\ = 200, m = 100). 
This behavior is described by a very simple approximated ex¬ 
pression valid in the limit C\ » e 2r » e 2 ', with not very 
large h\p, and when 6 , e —» 0 (corresponding to y, A <K k ), 

e 2r 

V ~ 2e _2r + (1 + hi + n 2 ) — (36) 

4C] 


which exhibits a minimum (hence corresponding to maximum 
entanglement) as a function of r at 


r ^ 



8 C 1 \ 

«1 + «2 + 1 / 


(37) 


given by v opt ~ . For values of r much larger or 

much smaller than this value, the resonators may not be entan¬ 
gled. When r is increased to very large values r » r opt , Q is 
reduced and the cooling dynamics becomes slow as compared 
to the standard mechanical dissipation, which takes place at 



0 50 100 150 200 

time (units of 1 /k) 

FIG. 2: Time evolution of the logarithmic negativity E N starting from 
an initial uncorrelated state with the optical mode fluctuations 8a in 
the vacuum state and each MR in its thermal state with mean phonon 
number: i) h\ = «2 = 0, Gi = 0.995G2 (black line); ii) h\ = 200, 
h 2 = 100, Gi = 0.918G2 (blue line); iii) h { = 1000, h 2 = 500, 
Gi = 0.82G 2 (green line); = 2000, n 2 = 1000, G f = 0.75G 2 (red 
line); the other parameters are y = 10 s' 1 , k = 10 5 s' 1 , G 2 = 10 5 s' 1 , 
A = 10 3 s' 1 . 
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rate ~ y (hj + 1), so that the correlations between the MRs 
cannot be efficiently generated. On the other hand, at small 
r <s r opt the Bogoliubov modes are essentially equal to the 
original modes, so that the cavity cools only the second res¬ 
onator, and also in this case mechanical entanglement can not 
be observed. Fig.[3]also shows that the simplified expression 
of Eq. ( [36] ) provides a simple but valid approximation for large 
C i and a very good estimate of the optimal value of the two¬ 
mode squeezing parameter r, i.e., of G 1 /G 2 , given by Eq. m 
The corresponding value of the logarithmic negativity is 


E n 


l ln [ Cl 
2 2 (1 + h\ + F 12 ) 


(38) 


and shows that once that the ratio G 1 /G 2 is optimized, the 
achievable stationary entanglement between the two MRs in¬ 
creases with increasing C\!{h\ + /I 2 ). 



r 


FIG. 3: E n at the steady state versus r for y = 10 s -1 , k = 10 5 s -1 , 
Gi = 10 5 s -1 , A = 0, implying a cooperativity C 1 = 2x 10 4 , and h\ = 
200, h 2 = 100. The full red line refers to the steady state solution 
of the QLE in Eqs. fBHTT), that is given by Eq. GD- the blue 
dashed line is evaluated with the approximated value of v reported in 
Eq. ( |34|), a nd the black dashed line corresponds to the approximation 
in Eq. (36). 

The above analysis of the stationary entanglement of the 
two MRs extends the results of Ref. 121 in various direc¬ 
tions. First of all, our model extends to the case of nonzero de¬ 
tuning A a model discussed in the Supplementary material of 
Ref. l3l . We see that a nonzero detuning has a limited effect 
of the dynamic of entanglement generation, providing only an 
effective increase of n^ 00 *, which however becomes negligible 
as soon as A <s k (see Eq. ©). Moreover, Ref. 01 provided 
an explicit expression for E N only for the case of negligible 
thermal occupancies and not too large values of r, while the 
present discussion applies for arbitrary values of r, h\ and h 2 . 


IV. DYNAMICAL EVOLUTION IN THE CASE OF EQUAL 
COUPLINGS 

In the special case of equal couplings G\ = G 2 = G, i.e., 
Q — 0, the Bogoliubov modes cannot be defined anymore and 
the description of the preceding Section cannot be applied. 
The dynamics is nonetheless interesting and still allows for 


the generation of appreciable entanglement between the two 
MRs, even though only at finite times and not in the station¬ 
ary state. We notice that in this special case, our scheme be¬ 
comes analogous to that of Ref. HI . that showed that two 
appropriately driven optical modes can be entangled with a 
pulsed scheme by their common interaction with a MR. More 
precisely, the QLE of Eqs. 0-0 are the same as those stud¬ 
ied in Ref. HU but now referred to two mechanical modes 
coupled to the same optical mode, i.e., with exchanged roles 
between optical and mechanical degrees of freedom. 

The physical mechanism at the basis of the generation of 
dynamical entanglement can be understood by looking at the 
Hamiltonian evolution of the system at equal couplings. Such 
mechanism essentially coincides with the one proposed for 
entangling internal states of trapped ions by Milburn ED and 
by Sprerisen and Mplmer H3l . and first applied to an optome¬ 
chanical setup by Kuzyk et al. M- In the present case, the 
common interaction with the bichromatically driven optical 
mode dynamically entangles the two MRs, and at special val¬ 
ues of the interaction time the optical mode is decoupled from 
the two MRs and mechanical entanglement can be strong. 

At equal couplings it is convenient to rewrite the effec¬ 
tive Hamiltonian after linearization of Eq. © in terms of 
mechanical and optical quadratures, using the expressions 
5bj = (xj + ipj)/ V2, j = 1,2, and 5a = (X + iY)/ V2. One 
gets 

H c „ = ^ (A 2 + Y 2 ) + hG V2 (x + X - p-9), (39) 

where jc± = (jti ± x 2 )/ V2, p ± = (p\ ± P 2 )/ V2 are linear com¬ 
binations of the two position and momentum operators of the 
two MRs. The Heisenberg evolution of these latter mechan¬ 
ical operators can be solved in a straightforward way, by ex¬ 
ploiting the fact that x+ and p are two commuting conserved 
observables. One gets (see also Refs. gaga ED) 


x+(t) = x+(0), p-(t) = p-( 0), (40) 

2 G 2 

x_(f) = iL(0) + — r (sinAr - At) p-(0) 

A 2 

2G 2 

+ ^(1 -cosAf)x + (0) (41) 

A 2 

gV 2 - gV 2 

-sin ArT(0) h—( 1 - cos Af)A(0), 

2 G 2 

P+(t) = p + ( 0) - (sin At - At) x+(0) 

9 

+ ^(l-cosAf)/3_(0) (42) 

GV2 - GV2 

-sin AfA(0)--— (1 - cos At) T(0). 


Relevant interaction times are those when the MR dynamics 
decouple from that of the optical cavity, and this occurs at 
t m = 2»i7t/A, m — 1,2,..., where 

2G 2 

x-(tm) = xJO) - 2mn-^-p-(0), (43) 


2 G 2 

p+(t) = p+( 0) + 2 mn-^-x+( 0 ). 


( 44 ) 












This map describes a stroboscopic evolution in which the 
two MRs become more and more entangled, because it cor¬ 
responds to the application of the unitary operator 


V. EFFECT OF THE COUNTER-ROTATING TERMS. 
STUDY OF THE EXACT DYNAMICS 


7 


U m = exp 
= exp 


2iTG~m 


A 2 


(*+ + P-) 


(45) 


2nG~m 


(6b\6b\ + Sb^pbn + 1 + 6 b\6b\ + 5b\6b2j 


This ideal behavior is significantly modified by the inclusion 
of damping and noise, especially the one associated with the 
cavity mode, which acts on the faster timescale 1 / k and seri¬ 
ously affects the cavity-mediated interaction between the two 
MRs, as soon as k becomes comparable to A. Mechanical en¬ 
tanglement is large for large G/A and we expect well distinct 
peaks for /Tv at interaction times t m , in the ideal parameter 
regime G » A » k. In the more realistic regime in which 
G, A and k are comparable, the peaks will be washed out, 
but we still expect an appreciable value for the mechanical 
entanglement for a large interval of interaction times. This 
is confirmed by the numerical solution of the time evolution 
associated with the QLE shown in Fig. |4j which refers to 
the parameter set y — 10 s' 1 , k = 10 5 s' 1 , G = 10 s s' 1 , 
hi = 200, Tin = 100, and to three different values of the de¬ 
tuning, A = 10 3 s' 1 (black dashed line), A = 10 4 s' 1 (red full 
line), and A = 10 5 s' 1 (blue full line). We see that an appre¬ 
ciable value of E n (even though smaller than the one achiev¬ 
able at the same h\ and Tin after the optimization of G 1 /G 2 of 
the previous Section) is reached for a large interval of interac¬ 
tion times t. Therefore even at equal couplings (and nonzero 
detuning) one can entangle the two resonators with a pulsed 
experiment. Mechanical entanglement instead vanishes in the 
stationary state. 



0 20 40 60 80 

time (units of 1 Ik) 

FIG. 4: Time evolution of E N in the case of equal couplings for the 
parameter set y = 10 s' 1 , k = 10 5 s' 1 , G = 10 s s' 1 , n, = 200, 
h 2 = 100, and for three different values of the detuning: A = 10 3 s' 1 
(black dashed line), A = 10 4 s' 1 (red full line), and A = 10 5 s' 1 (blue 
full line). 


The derivation of the effective linearized dynamics of Ap- 
pendix[B]suggests that the counter-rotating terms that we have 
neglected may play an important role when the mechanical 
frequencies are not too large with respect to the other param¬ 
eters (see also the comments in the supplementary material of 
Ref. |34| ). It is therefore interesting to study their effect by 
comparing the above predictions, both in the case of G 2 > G\ 
and in the case of equal couplings, to the solution of the exact 
QLE obtained without neglecting the various time-dependent 
terms. 

In Appendix [B] we describe the derivation of the effective 
linearized equations that we have studied in the preceding sec¬ 
tions and that is based on the elimination of fast rotating terms 
and on the expansion of the linearized coupling strength at 
lowest order in gj. Here we analyze the limit of validity of 
these approximations by solving numerically the system dy¬ 
namics with the inclusion of the non-resonant terms expanded 
at different orders in powers of gj. In Fig. [5] and [6] the red 
lines are evaluated without the non-resonant terms (i.e., the 
treatment of the preceding Sections), while the green and the 
blue ones take into account the full dynamics. In particular 
the green lines are computed by expanding the average fields 
a(t) and fl/t) (that have been introduced in Eq. ([2j) and dis¬ 
cussed in Appendix |Bj, at the lowest relevant order in pow¬ 
ers of g. while for the blue ones they have been expanded 
at sixth order in powers of g. Moreover, the green line re¬ 
sults are found considering only the steady state solution for 
a(t) and f3j(t), while the blue lines are computed taking into 
account their full dynamics (that includes also the transient 
regime before the steady state is reached) with initial condi¬ 
tion a( 0) = Pj( 0) = 0. 

In Fig. [5] we compare the time evolution of the entangle¬ 
ment evaluated with and without the time-dependent terms 
when Gn> G\. The parameters used in these plots are consis¬ 
tent with those used in Fig. [2] Specifically the three red curves 
in Figs. [5] (a), (b) and (c), that are barely visible because al¬ 
most entirely covered by the green curves, are equal to the 
three lowest curves in Fig. [2] We observe that the green and 
the red lines are always very close, meaning that the linearized 
RWA treatment is a very good approximation of the full dy¬ 
namics when a(t) and Pj(t) can be expanded at lowest order 
in g. Nevertheless, we note that if the mechanical frequencies 
are not large enough and higher order terms are taken into ac¬ 
count together with the full dynamics of a(t) and f3j(t). then 
the results can be significantly different as described by the 
blue curves. Specifically, the solid-blue lines are evaluated for 
sufficiently large values of the mechanical frequencies so that 
the condition in Eq. (|3]i is well fulfilled, and the effective lin¬ 
earized RWA dynamics recovers with significant accuracy the 
one determined with the inclusion of the non-resonant terms. 
The dashed-blue lines are instead evaluated for smaller fre¬ 
quencies. In this case it is evident that the non-resonant terms 
have a significant role in the system dynamics and that the 
lowest order expansion of the coefficients a(t ) and does 
not provide an accurate description. We note that according to 
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FIG. 5: Comparison of the time evolution of E N evaluated with and without the non-resonant terms, when G 2 > G\. The red lines are evaluated 
with the model described by Eqs. 0-0 that does not take into account the non-resonant terms. The green lines are evaluated with Eq. {B23}, 
which takes into account the non-resonant terms by considering the expansion of the steady state solutions of a{t) and pj(t), at first order in g as 
defined in Eq. (B22j. The blue lines are evaluated instead with Eq. l|B19f by considering the expansion for off) and fijit), calculated iteratively 
with Eqs. ( |B6fr . ( |B9[ )-( |B 1 2| ), up to sixth order in g; in particular these results take into account the full dynamics of the average fields a(t ) and 
pj(t), with initial condition a(0) = j8 ; (0) = 0, and not only the steady state as in the case of the green lines. The solid lines refer to a> 2 = 100k 
and aii = 50k, while the dashed lines refer to at 2 = 50k and oq = 25 k. The other parameters are G 2 = k, A = 0.01k, y = 10~ 4 k, and g = 10~ 4 k. 
Moreover in (a) G i = 0.918k hi = 200, h 2 = 100; in (b) Gi = 0.82k, hi = 1000, h 2 = 500; in (c) G\ = 0.75k, hi = 2000, h 2 = 1000. 





FIG. 6: Comparison of the time evolution of E N evaluated with and without the non-resonant terms, when Gi = G 2 . The insets report the 
corresponding average photon number (6a^ 6a V Red, green and blue lines are evaluated as in Fig.jsj The solid lines are found for A = 0.01k 
and the dashed lines for A = 5k. In (a) y = 0.03k, w, = 58k, a> 2 = 100k, h\ = 2 and h 2 = 1; in (b) y = 0.01k, uq = 58k, a> 2 = 100k, hi =2 and 
h 2 = 1; in (c) y = 0.001k, oq = 51k, u> 2 = 100k, hi = 20 and h 2 = 10. The other parameters are G\ = G 2 = k and g = 10~ 4 k. 


Eq. Q, in order to eliminate the fast rotating terms, the ratios 
cOj/Gj have to be much larger than one. Although the dashed- 
blue curves are evaluated for a ratio of roughly 25, 

which can be considered significantly large, we have found, 
indeed, that it is not enough for a faithful approximation of 
the system dynamics with the model discussed in the preced¬ 
ing sections. The conclusive analysis of these cases would, 
possibly, require a non-perturbative approach that is beyond 
the scope of the present work. A final remark is in order. We 
have verified that the discrepancy between the dashed-blue 
lines and the red ones is due to the combined effect of the 
higher order terms and of the transient initial dynamics of a(t) 
and [ij(t). Specifically, when we consider either the lowest or¬ 
der terms and the transient dynamics, or the higher order terms 
and only the steady state of a(t) and fij(t), the corresponding 
results for the entanglement dynamics are very similar to the 
red lines. 

In Fig. |6] we study the case of equal couplings G i = GV 
In this case solid and dashed lines differ in the values of the 
cavity detuning A. In general larger A (dashed lines) corre¬ 
sponds to smaller entanglement, and the results evaluated by 


including the counter-rotating terms tends to exhibit larger en¬ 
tanglement than the corresponding ones obtained without the 
non-resonant terms. The solid curves are found with smaller 
A. In this case red, green and blue lines are very close when 
the mechanical dissipation is sufficiently large as in Fig. [0(a). 
Larger discrepancies are found when the mechanical dissipa¬ 
tion is reduced as in Fig.[6](b) and (c), especially at relatively 
large time. We observe in fact that, while the red curves for 
the entanglement decay to zero at large time, the correspond¬ 
ing green and blue lines seem to approach a finite sizable 
value. As shown by the insets, when this different behaviour 
is observed, the average photon number in the cavity tends to 
diverge. This is a signature of the fact that the full dynam¬ 
ics including counter-rotating terms is actually unstable, even 
though the RWA dynamics without these terms is stable (see 
Eq. 0). We have confirmed the unstable nature of the time- 
dependent dynamics by calculating the Floquet exponents of 
the dynamical equations of the system. In fact, when a{t) and 
f5j(t) are considered in their steady state, one has a system 
of linear differential equations with periodic, time-dependent 
coefficients (see Appendix [0, and the Floquet theory can be 
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applied in this case lf52l ; we have verified that for the param¬ 
eters of Fig. [6] there is always at least one positive Floquet 
exponent, meaning that the system is unstable. This implies 
that, in general, the corresponding results are well-grounded 
only for relatively short time until the populations are not ex¬ 
ceedingly large. On the other hand, our results show that in a 
pulsed experiment with the parameters of Fig. [6j these insta¬ 
bilities do not constitute a serious hindrance to the creation of 
significant entanglement at finite times. 

Therefore, when the mechanical frequencies are sufficiently 
large (toj > 10 2 k) (and, limited only to the case of equal cou¬ 
plings, when also mechanical damping is not too small), the 
effective linearized RWA dynamics obtained by neglecting the 
counter-rotating terms approximates with very good accuracy 
the full system dynamics. 


VI. STRATEGIES FOR THE EXPERIMENTAL 
DETECTION OF MECHANICAL ENTANGLEMENT 


We finally discuss how to detect the generated mechani¬ 
cal entanglement between the two MRs at different frequen¬ 
cies. The present entanglement describes EPR-like correla¬ 
tions between the quadratures of the two MRs and therefore 
we need to perform homodyne-like detection of these quadra¬ 
tures. In the linearized regime we are considering the state 
of the two MRs is a Gaussian CV state, which is fully char¬ 
acterized by the matrix of all second-order correlations be¬ 
tween the mechanical quadratures. Therefore from the mea¬ 
surement of these correlations one can extract the logarithmic 
negativity E^. One does not typically have direct access to 
the mechanical quadratures, but one can exploit the currently 
available possibility to perform low-noise and highly efficient 
homodyne detection of optical and microwave fields, and im¬ 
plement an efficient transfer of the mechanical phase-space 
quadratures onto the optical/microwave field. 

As suggested in Ref. f53l and then implemented in the elec¬ 
tromechanical entanglement experiment of Ref. a, the mo¬ 
tional quadratures of a MR can be read by homodyning the 
output of an additional “probe” cavity mode. In particular, if 
the readout cavity mode is driven by a much weaker laser so 
that its back-action on the mechanical mode can be neglected, 
and resonant with the first red sideband of the mode, i.e., with 
a detuning A p . — oj j, j = 1,2, the probe mode adiabatically 
follows the MR dynamics, and the output of the readout cav¬ 
ity a°j l " is given by (see Fig. |TJ) |53l 

G p . 

a'f = i -p 6bj + a 1 ", j = 1,2, (46) 

1 \[k ' 1 

with G P j the very small optomechanical coupling with the 
probe mode. Therefore using a probe mode for each MR, 
changing the phases of the corresponding local oscillator, and 
measuring the correlations between the probe mode outputs, 
one can then detect all the entries of the correlation matrix and 
from them numerically extract the logarithmic negativity E^. 


A. Concluding remarks 


We have studied in detail a general scheme for the genera¬ 
tion of large and robust CV entanglement between two MRs 
with different frequencies through their coupling with a single, 
bichromatically driven cavity mode. The scheme extends and 
generalizes in various directions similar schemes exploiting 
driven cavity modes 04l l36i [38] ;44) for entangling two MRs 
or two cavity modes. The scheme is able to generate a remark¬ 
ably large entanglement between two macroscopic oscillators 
in the stationary state, i.e., with virtually infinite lifetime, and 
it is quite robust because one can achieve appreciably large 
CV entanglement even with thermal occupancies of the order 
of 10 3 . The scheme is particularly efficient in the limit where 
counter-rotating terms due to the bichromatic driving of the 
cavity mode are negligible, and we have verified with a care¬ 
ful numerical analysis that this is well justified when the two 
mechanical frequencies are sufficiently large ajj > 1 0 2 k. 
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Appendix A: Normal modes and Hamiltonian dynamics 


It is straightforward to see that the diagonal form of the 


interaction Hamiltonian of Eq. (10 1 is 


H e fj = HAo0\0i + hAia\&i + hA2& 2 &2, (Al) 


where 


a i = cos 602 + sin 66a, 
&2 = cos 66ci - sin 602, 


(A2) 

(A3) 


define the other two normal modes together with the dark 
mode 0\, introduced in Eq. 0 , with 6 defined by the condi¬ 
tion tan 2 6 = -2Q/A, while the eigenvalues are given by do = 
0, Ai = (A - A) /2, A 2 = (A + A) /2, with A = 0A 2 + AQ-. 

The normal modes allows to understand the dynamics in 
the absence of optical and mechanical damping processes. 
In fact, from Eq. ( |A1| > one can easily derive the Heisenberg 
evolution of the mechanical bosonic operators. By invert¬ 


ing Eqs. (11 1 -( 12 1 one has, 6b\(t) = coshr/?i(f) - sinh r01(t). 


6b 2 (t ) = cosh r0 2 (t) - sinh r0 J .(t) and using aj{t) = e ut a0 0), 
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j = 1,2, and fii(t) = /?i(0), one gets 
5b\(t) = coshryg^O) 

- sinh r exp 


(A4) 


.A t 

At 

At 


cos —— 
2 

i cos 26 sin — 
2 


m 


-i sinh r sin 20 exp 


6 b 2 (t) = - sinhry§|(0) 
At 


At 


At 

sin — 6 a( 0) 


(A5) 


+ cosh r exp 


+i cosh r sin 26 exp 


At At 

cos-h i cos 26 sin — 

2 2 


> 02 ( 0 ) 


At 


At 

sin — 6 a( 0 ). 


We now look for special time instants at which the two 
mechanical modes can be strongly entangled. A necessary 
condition for such dynamical entanglement is that at these 
times, the cavity mode must be decoupled from the mechan¬ 
ical modes and Eqs. (A4 )-( |A5] > show that it occurs when 

sinAf/2 = 0, i.e., t p = 2np/ sjA 2 + AQ 2 , p = 1,2,_ At 

these time instants one has 


6 b\{t p ) = [cosh 2 r - sinh 2 r] 5Z?i(0) (A6) 

+ sinh/-coshr(l - e l<f>p ) 6 b\( 0), 

6 P 2 (tp) = [y^ cosh 2 r - sinh 2 r] 6 b\( 0) (A7) 

- sinhrcoshrfl - e'M<5bi(0), 

where = np(l + A/AV In particular, if = -1 one gets 

6 b\(t p ) = cosh2r<5Z?i(0) + sinh2r<J^J,(0), (A8) 

5b\(t p ) = - cosh 2r6b\(0) - sinh2r<5£>i(0), (A9) 


i.e., the state of the two MRs at time t p is the result of the 
application of the two-mode squeezing operator with squeez¬ 
ing parameter 2r, S (2 r) (see Eq. (14 1 ) to their initial state. In 
the usual case of an initial thermal state for the two MRs with 
mean thermal phonon numbers hj, the state at time t p is there¬ 
fore a two-mode squeezed thermal state ESI (see Eq. ©), 
with logarithmic negativity Il47li48l 


Efj(tp) = --In [« 2 + (n+ + l) 2 cosh8r (A10) 

- yl(h+ + l) 4 sinh 2 8r + 4 n 2 _(n + + l) 2 cosh 2 4r , 


where h ± = n\ + hi. For the relevant case of not too small 
values of the squeezing parameter r, E N can be well approx¬ 
imated with its value at equal mean thermal phonon number 
h- = 0 , 


EN(tp) — 4r - In [n + + 1], (All) 

showing that at this interaction time, the entanglement be¬ 
tween the MR can be very large, even if starting from a rela¬ 
tively hot state, by properly tuning the ratio G 2 /G\, i.e., the in¬ 
tensity of the two tones. This large mechanical entanglement 


is achieved when the condition = -1 is also satisfied for 
a given integer p. This is obtained for any odd p when A = 0, 
or by properly adjusting the value of Q for a given A 4- 0, i.e., 
if 


Q 1 = A 2 ^—^^ d odd, 0 < d <2p, d + p. (A12) 
' 4 (p - dY 

This dynamical scheme for the generation of continuous vari¬ 
able mechanical entanglement is similar to the Bogoliubov 
scheme proposed in Ref. B31 for entangling two optical cavity 
modes. It is extremely hard however to use it for entangling 
two mechanical modes as in the present case, because the cav¬ 
ity decay rate is comparable to Q and A in typical situations, 
thereby strongly affecting the ideal Hamiltonian dynamics de¬ 
scribed here. 


Appendix B: Linearization of the optomechanical dynamics 
with two-frequency drives 

The system dynamics is described by the following QLE 
a = - [k + i (Ao + co-)] a - i [fije^'"*' + £ , 2 e'“ +( ] + V2 ~Ka' n 

-i [gi (b\ + b\) + g 2 (b 2 + y)l 3. 

bj = - (y + icojJ hj - igjtfci + y/Yjb'f j = 1,2, 

(Bl) 

where, here, differently from the description used in Sec. [H] 
we are representing the cavity field in a reference frame rotat¬ 
ing at the frequency a>o -(to 2 - to\)/2, and we have introduced 
the frequencies 


(j)2 ± (j) \ 


The other parameters and operators are defined in the main 
text. 

If we perform a time dependent displacement, for both cav¬ 
ity and mechanical degrees of freedom, of the form 

d(t) = 6 a(t) + a(t), 

bj(t) - Sbj(t) + Pj(t), (B2) 


the QLE reduce to the form 

6 a = -{K + i (Aq + w_) + 2/giRe [/?i(0] + 2 ig 2 Rc \{3 2 (t)]\6a 


—iA(t) + V2 ~Ka m 


-ia(t) [gi (shi + 6b \) + g 2 ( 6b 2 + h^)] 


-i [gi (bhi +6b\) + g 2 (Sb 2 + 6b\^ 6a, 

(B3) 

6b j = - (y + iajjJ 6b j - iBj(t) + y/f]bj 


-igj ]a(t)6a ! + a(f)*daj - igj6a^6a. 

(B4) 
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where the new driving terms, A(t ) and Bj(t) read 
A(t) = + E 2 e lb>+t } - - [ a : + i(A 0 + a>-)]a(t) 


-2 i (giRe |j3i(0l + g 2 Re \[h(t)]) a(t), 

d/3j(t ) /y,- \ , 

B/t) = - (^ + iw ; )j8 y -(0 - igj |or(OI“ • 


When gj and go are sufficiently small and we chose a(t) and 
Pj(t) such that A(t) = 0 and Bj(t) = 0, then the non-linear 


terms, i.e. the last terms in the two equations (B3 i and (B4i, 
can be neglected. The equations A(t) = 0 and Bj(t) = 0 define 
a set of non-linear differential equations with periodic driving 
for the parameters off) and [3jit). The solution can be eval¬ 
uated perturbatively in the small parameters g\ and go 1461 . 
Here we assume gi - g2 = 8 and we observe that the solu¬ 
tions for a(t) and ftjit), with initial condition a(0) = fijiO) = 0, 
contain, respectively, only even and odd powers of g. 


ait) = 'Yj S P a (p \t), 


P = 0 
p even 


Pjit) = 2 s p P { f\t). 


(B6) 


p = i 

p odd 


The equations for each component of these expansions can be 
written in the form 


a (p \t) = -za yp> {t) + '3£ > {i), 
Pfit) - + S (p \t), 


( P)(A 4- E'O’V 


where 


z = k + i (A 0 + w_), 

7j 4. ■ 

Wj = ~+lCOj, 

and the driving terms are defined recursively as 

p -1 


(B7) 


(B8) 


S Tit) = -2 i^a^iO 


9=0 


x (Re \p ( r H \t)\ + Re \p ip ~ q ~ V \t )\), 


ajfV) = -i^a^it) a^-\f) , 

q= 0 

with the initial condition 

S^ 0) (f) = E 1 e~ i0>+t + E 2 e ib>+t , 

aj®(0 = o. 


(B9) 


(BIO) 


In particular they can always be rewritten as sums of exponen¬ 
tial functions of the form 


^(f) = 

n 


with x P ' n \ x) P ' n) i ^a' n) and time-independent complex 
coefficients, whose specific form can be computed iteratively. 
Moreover, the expression for a (p \t) and fi (p> (t) are found inte¬ 
grating Eq. and are given by 


(B5) 


a 


w « - zd^ K'-i- 


z+<r; 


/*?’(') = Z K”'" e ""'1 ■ (B12) 

■P 


n Wj + C 


We note that all the coefficients an d Cp have non- 
positive real parts. Re [4of’ n) ]. Re < 0, thus the large¬ 

time solutions cr (p) (f — > oo) = a (p \t) and [3 {p \t — > oo) = pff^it) 
are found from Eq. ( |B12[ ) by keeping only the terms for which 
i^ p ' n> and ^ p ' n) are purely imaginary, that can be shown to be 
equal to i n a>+ with n odd and even integer respectively. In 
particular a (p \t) and p {p ] t it) are periodic functions (with pe¬ 
riod 2 nlu> + ) which contains frequency components that are, 
respectively, odd and even multiples of a> + . 


p +1 

oJ'W = Z 

n = —p — 1 
n odd 


x, 


(p,n) 


Z + i n u> + 


J ncj+t 




p +1 

z 


X P 

i Wj + ina)+ 


J ncj+t 


(B13) 


X (p ' n \t) are the coefficients that correspond to those particu¬ 
lar parameters £ pX) and i^' n) that are imaginary. 

1. Resonant and non-resonant terms 

The QLE, in the interaction picture with respect to the 
Hamiltonian Hq = h(u>-cPa + <x>ib\b\ + , reduce to 


where z and vv ; are defined in Eq. (|B8b, and x^ p ’ n> it) and 


6 a = -(k + ;Aq) 6 a + V2 ~Ka m 

+2 i (giRe \Piit)] + g 2 Re \(h{t)]) 6a 
-ia(t)e lbJ -' [gi (b'De 
+g 2 ( 6 b 2 e- WJ2 ' + 6 bU^-')\ , 

Vj 
2 


(B14) 


~ imt + 6b\e mt ) 


6bj = -y 6bj+ ^Ifjbj - igj<z u ' lj ‘ \e laj -'ait)6a f + H.c. | . 


Before proceeding, we note that we can include the DC com¬ 
ponent of fij(t) into the cavity detuning, hence we introduce 


A = A 0 + 2 Xj= i,2 8/Re 

introduced in Eq. & 


yW» 

Y nP X J— 

ZjpS w , 


, according to the notation 


y(P -°) 

z*'— =i>f 

W : J 


(Bll) 


(B15) 
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Moreover we can isolate the resonant terms of the QLE, 
namely the terms with time-independent coefficients, by con¬ 
sidering the lowest order frequency components of a(t), i.e.. 


a± 



X. 


(p,± i ) 


Z+ICO+ 


(B16) 


min {o»i, W 2 , |wi - a» 2 |} ■ In particular this condition is true 
when it is valid for the lowest order term in the expansion in 
power of g. In details, the non resonant terms can be neglected 
when 


g |a ± |, k «: min{tui,w 2 ,|a>i - a> 2 \} . (B21) 


corresponding to the frequencies ±a>+, and defining 


Pjit) 

- Pjit)~Pf, 

(B17) 

a+i t) 
a-it) 

= ait) - e'" +f a+, 

= ait) - e" ! " +, a_ . 

(B18) 


Thereby we find 

6a = -(k + ;A) 6a - ia-g\5b\ - ia+g26b2 
+ 02~Ka m + F a (t ), 

6b\ = -~6b\ - igia-6aP + \[y{b™ + F bl it), 

6b 2 = -^6b 2 -ig 2 a + *6a+ rfFb' 2 n + F b2 (t), (B19) 

where F a (t), F bl (t ) and F bl it) account for the terms with time- 
dependent coefficients and are given by 

F a (t) = +2/(giRe + g 2 Re [/j 2 (f)]) 

-ia-(t)e io> -‘gidb\e io> ' 1 2 3 4 5 - ia + (t)e m - t g26b 2 e- im \ 
F bl (t ) = -/gie'" 1 ' + e~ ,w - t a*it)6a^ , 

F bl (t) = -ig 2 e WJ2> a(t)6a f + e~ l ' J> - , a + (t)‘ 5a\ . (B20) 


In particular we can introduce the linearized coupling st rength 
G ( [ ot) = ga _ and G+ ol) = ga+, with a ± defined in Eq. (BI 61 . 
The expressions introduced in Eq. 0 correspond to the ex- 
pansion of these parameters at zeroth order in g (see also 
Eq. ( |B22| ). 

We are interested in the regime in which the terms in 
Eqs. (B20) with time-dependent coefficients are negligible. 
They can be neglected when gj \a st (t)\,gj LS;, s r(f)j, k «: 


When this condition is fulfilled the parameters a(t) and fijit) 
can be safely expanded at the lowest order in g. Specifically 
they can be approximated as 


a — ^ 
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-iE 2 
z + icj+ 


a s ,(t) ^ og. it) = as ItJ+ ' + a+e'" +f . 


0 


iDC 


l Sj r * 
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• a+a* + ], 




gJ0%(O ~0f 


= ~ l §.i 


wj - 2ia>+ 


-2 i cj+t _ 


a + a_ 


(B22) 


wj + 2 ia>+ 


a 2i 01+ t 


Moreover the parameters a±(t ) defined in Eq. ( B1 8| ) are zero. 
Using these expressions the QLE in Eq. ( B1 9| ) can be rewritten 
as 


6a = - (k + /A) 6a - iGi6b\ - iG26b 2 
+ V2 ~Ka in + F a (t ), 

6 b\ = -2^-6bi - iG\6a^ + yfyib'i + F bl (t), 

6 b 2 = ~6b 2 - iG\6a + Vtt&T + F h (t), (B23) 
with G\ and G 2 defined in Eq. (j7]) and 

Fait) - + 2 i(^Rep 1 >)]+^Re[^( 0 ])tfa 
-i^-'a^it) [g x 6b^ t + g 2 6b\e “«), 

Fb.it) =; -igi^'+^affifida, 

F bl it ) - -ig 2 e i(uM a^it)6a^. (B24) 


When the time-dependent coefficients are neglected these 
equations reduce to Eqs. 0-0- 
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